Direct numerical simulation of dispersed particles in a compressible fluid 
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We present a direct numerical simulation method for investigating the dynamics of dispersed 
particles in a compressible solvent fluid. The validity of the simulation is examined by calculating 
the velocity relaxation of an impulsively forced spherical particle with a known analytical solution. 
The simulation also gives information about the fluid motion, which provides some insight into the 
particle motion. Fluctuations are also introduced by random stress, and the validity of this case is 
examined by comparing the calculation results with the fluctuation-dissipation theorem. 



I. INTRODUCTION 

Particle dispersions have various unique properties, 
and an understanding of these properties is important 
in many fields of science and engineering. These proper- 
ties originate from the dynamics of particles, which are 
extremely complicated because of the hydrodynamic in- 
teractions among particles mediated by the motion of 
the surrounding fluid. Therefore, several numerical ap- 
proaches have been formulated to investigate the dynam- 
ics of such dispersions. As one of these approaches, a 
direct numerical simulation has been developed, wherein 
the hydrodynamic interactions arc directly computed by 
simultaneously solving for the motion of the fluid and the 
motion of the particle. In recent years, we have devel- 
oped an efficient direct numerical simulation scheme for 
dispersions, which is called the smoothed profile method 
(SPM) [il, 0] . In this scheme, the sharp interface between 
the fluid and the particles is replaced by a smoothed in- 
terface using a continuous profile function. The Navier- 
Stokes equations are solved for the fiuid motion on a fixed 
square grid, and Newton's and Euler's equations of mo- 
tion for the particles are solved simultaneously while con- 
sidering the momentum exchange between the fiuid and 
the particles. 

The hydrodynamic interactions are transmitted in two 
ways: via viscous momentum diffusion and via sound 
propagation. The time scale of viscous diffusion over the 
particle size is t^ — a? /v, and that of sound propaga- 
tion is Tc = a/c, where a is the particle radius, v is the 
kinematic viscosity, and c is the speed of sound in the 
fiuid. Here, we will define the compressibility factor as 
the ratio of the two time scales: 



Ti, ac 



(1) 



This factor represents the degree of infiuence of com- 
pressibility in the dynamics of the dispersed particles. 
According to Eq. ([T]), the compressibility becomes in- 
creasingly important as the particle size decreases. In- 
deed, the compressibility of a system has been consid- 
ered an important factor in molecular scale dynamics 
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of monatomic liquids studied using hydrodynamic the- 
ories 0, Q. On the other hand, in studies of particle 
dispersions, compressibility has rarely been considered, 
and an incompressible host fiuid is generally assumed. 
In the case of a dispersed particle of radius a = 100 nm 
in water, for instance, the compressibility factor is evalu- 
ated from i^ = 1.0 X 10~^m^/s and c = 1.5 x lO^m/s to be 
£ = 6.7x10^"^. In many cases, researchers are interested 
in phenomena progressing over the time scale of viscous 
diffusion or even longer time scales, such as those relating 
to shear properties, electrophoresis, and sedimentation. 
Therefore, the assumption of an incompressible host fiuid 
is valid, and most direct numerical simulation methods, 
including SPM, have been developed on the premise of 
incompressible fiuids. However, when we investigate phe- 
nomena associated with sound propagation, such as sonic 
agglomeration, acoustic spectroscopy, and electroacous- 
tic measurements, the consideration of compressibility is 
required. 

In the present study, we extend the SPM to compress- 
ible fiuids. Some aspects of the dynamics of a single 
particle in a compressible fiuid have been theoretically 
analyzed ^Wt, and we compare the simulation results 
obtained herein with analytical solutions to determine 
the accuracy of the simulation. In particular, we con- 
sider the velocity relaxation of a spherical particle after 
an impulsive force is added. The numerical simulation 
also gives information regarding fluid motion for which 
the analytical solution is unknown, and the dynamics of 
the particle can be investigated from the viewpoint of the 
fiuid dynamics. In addition, we also consider a system 
with thermal fiuctuations by introducing random stress, 
and the velocity autocorrelation function is compared 
with the analytical solution according to the fiuctuation- 
dissipation theorem. 



II. SIMULATION METHOD 

A. Equations 

In the SPM, the particle-fiuid boundary is replaced by 
a continuous interface. For this purpose, the smoothed 
profile function </>(»", i) G [1,0] is introduced. This func- 
tion represents the boundary between the fluid and par- 



tide regions: (/> = 1 at the particle domain, and (j) — 
at the fluid domain. The two regions are smoothly con- 
nected through thin interfacial regions with a thickness 
characterized by ^. The mathematical expression of cf) 
can be found in a previous paper [Ij. The total velocity 
field is defined as 



(1 - (j))vf + 4)1 



(2) 



where Vf{r,t) is the fluid velocity field, and Vp{r,t) is 
the particle velocity field constructed from rigid motions 
of the particles [l|, |2|- When a compressible host fluid 
is considered, the fluid mass density is altered, and we 
define the mass density field as 



P/ = (1 - 0)P- 



(3) 



The auxiliary mass density field p{r, t) is defined over the 
entire domain. However, the physical fluid mass density 
Pf{r, t) must be zero within the particle domain, and this 
requirement is satisfied by multiplying by (1 — 0). 

The equations governing the dynamics of the disper- 
sion system are given as hydrodynamic equations with 
the addition of a body force. The hydrodynamic equa- 
tions consist of three conservation laws concerning mass, 
momentum, and energy. The conservation equations of 
mass and momentum are given as 



dp 

dt 



V • m = 0, 



dm _ , , „ 

— + V • (mv) = V • cr + p(j)fp, 



(4) 



(5) 



where m{r,t) — p{r,t)v{r,t) is the momentum density 
field. We consider a compressible Newtonian fluid, and 
the stress tensor is given by 



-pi + r][Vv + (Vi;)^ 



Vv 



ry (V • v)I,{6) 



where p{r, t) is the pressure, rj is the shear viscosity, and 
77„ is the bulk viscosity. A body force p4>fp is added so 
that the rigidity of the particles is satisfled. Additionally, 
we assume a barotropic fluid described by p = p(p), and 
the pressure gradient is proportional to that of the mass 
density: 



Vp = (?Vp. 



(7) 



Equations (H])-© are closed for the variables p, m, and 
p; therefore, energy conservation does not need to be 
considered. 

The motion of the dispersed particles is governed by 
Newton's and Euler's equations of motion: 



M.^V.^F- 



^R^ - F„ 



(8) 



where i?i, Vi, and il^ are the position, translational ve- 
locity, and rotational velocity of the i-th particle, respec- 
tively. The particle has a mass Mi and a moment of iner- 
tia li. The hydrodynamic force F^ and torque N^ are 
exerted on the particle by the fluid, and the force F^ is 
exerted through direct interactions among the particles. 
The effect of thermal fluctuations on the dynamics of 
particles is important when the particle size is on the 
order of a micrometer or smaller. Fluctuations are intro- 
duced through a random stress tensor s, which is added 
in Eq. ([5]). The random stress is a stochastic variable 
satisfying the fluctuation-dissipation relation [8] : 

(szj(r,t)sfe,(r',t')) = 

2kBTTj,,kiS{r' ~ r)5{t' - t), (10) 

where fc^ is the Boltzmann constant, T is the tempera- 
ture, and 



Vijki = vi^ikSji + SaSjk) + ( Vv - i^V 1 S.jSki- (11) 

Brownian motion of the dispersed particles is induced by 
random stress acting on the fluid. Thermal fluctuations 
can be introduced using the Langevin approach, where 
random forces are exerted on the particles [9!]. However, 
this approach does not accurately represent the short- 
time dynamics of the system because the time autocorre- 
lation of the hydrodynamic force acting on the particles 
is neglected. Therefore, the fluctuating hydrodynamics 
approach is more appropriate for investigating dynamics 
over the time scale of sound propagation. 



B. Simulation procedure 

Here, the time-discretized evolution of the equations is 
derived. The time step i„ represents the n-th discretized 
time, and the time step change from i„ to t„+i = tn + h 
will be considered. The time evolution of the system is 
determined through the following steps: 

(i) The mass and momentum density changes associ- 
ated with sound propagation are calculated as 



n+l n 

P = P 



t„+h 



dtV ■ m, 



tn 



m* =m^ - c^ 



tn+h 



AtVp. 



(12) 
(13) 



When we assume a periodic boundary condition 
and use the Fourier spectral method, a semi- 
implicit scheme becomes feasible [l3|- This situ- 
ation eases the restriction on the time increment 
for a small compressibility factor e. 

(ii) The time evolution of the advection and viscous 
diffusion terms are calculated as 



dt 



n, = AT. 



H 



(9) 



tn+h 



diV • {t — mv) 



(14) 



where r is the dissipative stress given by <t = —pl- 



T. 



(iii) In concert with the advection of the particle do- 
main, the position of each dispersed particle evolves 
as 



R 



>n-\-l jyn 



tn-\-h 



dm. 



(15) 



(iv) The hydrodynaniic force and torque are derived by 
considering the conservation of momentum. The 
time- integrated hydrodynamic force and torque are 
computed as 



rtn+h 



I diJ;^ = y drC+^m" - p"+i<), (16) 



t„+h 



dtN, 



H_ 



Jdrlir^Rr') X ^r\rn** ^ p-+'v;)]. (17) 



With these and other forces acting on the particles, 
the translational and rotational velocity of each dis- 
persed particle evolve as 






dtiFf+Ff), (18) 



III. SIMULATION RESULTS 

Numerical simulations are performed for a three- 
dimensional box with periodic boundary conditions. The 
space is divided by meshes of length A, which is the unit 
length. The units of the other physical quantities are de- 
fined by combining t] — 1 and po = 1 with A, where po is 
the fluid mass density at equilibrium. The system size is 
Lx X Ly X Lz — 128 X 128 x 128. The other parameters 
are set as a = 4, ^ = 1, pp = 1, rjy = 0, and h = 0.01, 
where pp is the particle mass density. We performed sim- 
ulations of the dynamics of an isolated spherical particle 
in a fluid in two situations. First, we investigated the 
relaxation response of a particle with an impulsive force. 
Second, we consider the velocity autocorrelation function 
of a particle with thermal fluctuations. 

Because the input particle radius a = 4 is not neces- 
sarily equal to the effective hydrodynamic radius a* of 
the particle represented by the smoothed profile func- 
tion (/) with a fuzzy interface of thickness ^, we calcu- 
lated a* from the drag force acting on the spherical par- 
ticle moving at velocity V, which is analytically given by 
Fd — 6TTr]aK{(j))V, where K{cl)) represents the effect of 
the periodic boundary condition depending on the vol- 
ume fraction of the particle ip lUl ■ The effective radius 
was evaluated as a* — 3.87 in the present simulations, 
and therefore, the momentum diffusion time is estimated 



to be Ti, 



/v in this case. When we compare the 



present simulation results with the analytical solutions, 
the corrected particle density p* = {a*/a)^pp and the 
compressibility factor e* = {a/a*)e of the analytical so- 
lutions are employed. 



Relaxation 



f2"+i =n" + /ri 



t„+h 



dtN, 



H 



tn 



(19) 



(v) The updated velocity of the particle region is im- 
posed on the velocity field as the body force pcfifp. 



m ^ = m 



tn + h 



dtp(f>fi 



P' 



(20) 



Let us consider a single spherical particle in a fluid at 
rest. We will investigate the relaxation of the particle 
velocity after exerting an impulsive force at the center of 
the particle. The impulsive force is assumed to be small; 
in other words, we will consider the motion of the par- 
ticle and the fluid for low Reynolds and Mach numbers. 
We set the impulsive force to achieve an initial particle 
Reynolds number of RCp = 0.01. Here, we introduce the 
velocity relaxation function as the normalized velocity 
change of the particle as 



tn+h 



dtpcf>fp = 0"+i(p"+it;;+i - m**). (21) 



In the case of an incompressible fluid, the pressure is 
spontaneously determined by the solenoidal condition of 
the velocity field. On the other hand, in the present case, 
the pressure or mass density varies independently of the 
velocity field. 



Vit) = -7(t), 



(22) 



where P is the impulsive force added at t = 0. The an- 
alytical form of the relaxation function is obtained from 
a Stokes approximation (see the Appendix). 

The simulations are performed with compressibility 
factors of £ = 0.1, 0.6, 1.0, and 1.5. The condition 
e > 1 implies that the sound propagation occurs more 
slowly than the viscous diffusion, according to the defini- 
tion given in Eq. ([T]). The simulation results agree quite 
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FIG. 1. (Color online) The velocity relaxation function for 
various compressibility factors: (a) e = 0.1, 0.6, (b) 1.0, and 
1.5. The bold lines illustrate the analytic solutions for various 
compressibility factors, and the bold dashed double-dotted 
line shows the analytic solution for an incompressible fluid. 



well with the analytical solutions for all of the compress- 
ibility factors e, as shown in Fig. [T] The oscillation of 
t/''':^ > 2 at e = 0.1 arises from the periodic boundary 
conditions. In this case, the sound pulse arrives at the 
end of the system at t/Ti, = 1.6; afterward, the sound 
pulse returns and affects the particle motion. 

The velocity relaxation results highlight some remark- 
able properties of dynamics in a compressible fluid. In 
the case of an incompressible fluid, part of the parti- 
cle momentum is instantly carried away by the propa- 
gation of the infinite-speed sound wave, and the parti- 
cle moves as if its mass were M* = M + Mf/2, where 
Mf = 4:Tra^po/3 is the mass of the displaced fluid [ia |. 
Therefore, the velocity relaxation function for an incom- 



pressible fluid at the initial time is 7(+0) = M/M*; af- 
terward, the particle velocity gradually decreases because 
of momentum diffusion in the fluid caused by its viscos- 
ity. Eventually, this decay obeys the power law t~'^/^ 
in the long-time region. On the other hand, in the case 
of a compressible fluid, the velocity relaxation function 
indicates that the relaxation due to sound propagation 
occurs in a finite time interval. For the compressibility 
factor e < 1, the two relaxation processes are almost sep- 
arate, and the relaxation function coincides with that of 
an incompressible fluid after the relaxation due to the 
sound propagation occurs. With an increase in the com- 
pressibility factor, nonmonotonic behavior is observed in 
the relaxation function, and finally, inversion of the par- 
ticle velocity is observed, as shown in the investigation 
of the analytical solution 7] • 

The density deviation around the particle is shown in 
Fig. [21 for which an analytical form has not been ob- 
tained. For a small compressibility e — 0.1, the sound 
wave pulse expands from the particle very quickly, which 
corresponds to a rapid decrease in the particle momen- 
tum due to sound propagation. On the other hand, for 
a large compressibility s — 1.5, the sound pulse does not 
spread. The pulse remains within the vicinity of the par- 
ticle and gradually decays through viscous diffusion. The 
continuation of high fluid density in front of the particle 
is expected to cause backtracking motion. 

The velocity field around the particle is also affected 
by the compressibility. The velocity field can be for- 
mally decomposed into two components: an incompress- 
ible (solenoidal) component v^ and a compressible com- 
ponent D*-". The former is a solenoidal vector field when 
'V ■ v^ — Q, and the latter is an irrotational vector field 
when V X t>'^ = 0. Therefore, the divergence of the 
velocity field V ■ v gives information regarding the com- 
pressible component, and the rotation of the velocity field 
V X t> gives information pertaining to the incompressible 
component. The rotation and divergence of the velocity 
field are shown in Figs. [3] and SI respectively. The ve- 
locity field is also depicted in these figures. An obvious 
vortex ring is observed around the particle for a small 
compressibility s = 0.1, while the vortex ring is not clear 
for a large compressibility e = 1.5. The rotation corre- 
sponds to the intensity of the vorticity, and the map in 
Fig. [3] represents a pair of vortex rings that is inherent to 
an incompressible fluid [l3| . The rotation remains nearly 
constant, regardless of the compressibility, which indi- 
cates that the incompressible component is not affected 
by the compressibility. The vector fields v^ and v'-^ in- 
fluence each other only through the nonlinear terms in 
Eqs. (|4]) and ([5]); however, in the present simulations, a 
low Reynolds number flow is assumed, and the additiv- 
ity of the incompressible and compressible components 
is almost valid. Therefore, the effect of compressibility 
is observed only from the divergence of the velocity field 
depicted in Fig. S) The region of positive divergence rep- 
resents the source of the flow, and the negative region 
represents the sink of the flow. From the equation of 
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FIG. 2. (Color online) Time evolution of the fluid density deviation around the particle. The compressibility factors are (a) 
e = 0.1 (upper row) and (b) e = 1.5 (lower row). The color scale (gray scale) represents negative (darker) to positive (lighter) 
density deviation. The black circle represents the particle. The direction of the initial particle velocity is right in these pictures. 



continuity Eq. (|3]), the divergence is equal to the reverse 
sign of the time change rate of the mass density: 



V • tJ = In p. 

Dt ' 



(23) 



According to this relation, the source corresponds to the 
density decrease and the sink corresponds to the den- 
sity increase. Therefore, the source and the sink move 
according to the sound propagation. The propagation 
speed decreases with increasing compressibility, which 
corresponds to the difference in Fig. |4] related to the 
compressibility. The pattern of the total velocity field is 
described as the superposition of the vortex convection 
of the incompressible component v^ and the source-sink 
flow of the compressible component v'^ . The compress- 
ibility factor governs the relative time evolution of each 
component to produce various flow patterns. 



B. Fluctuation 

Thermal fluctuations are introduced through random 
stress in the host fluid. Computationally, a random stress 
term satisfying Eq. (|10p is added to the stress tensor 
given by Eq. ([6]). We consider a single particle in the 
fluid, which moves randomly as a result of thermal fluc- 
tuations. From the fluctuation-dissipation theorem, the 
velocity autocorrelation function of the particle is related 
to the relaxation function as 

l{t)^^{V{Q)-V{t)). (24) 

The accuracy of the fluctuating system can be confirmed 
by the validity of this relation. In the numerical pro- 



cedure for solving a stochastic differential equation, the 
numerical error can be larger than that in an ordinary 
differential equation due to the truncation error in the 
time integration of the random noise term \\4\ . This er- 
ror is decreased with the decrease in the time increment, 
and we set the smaller time increment than that in the 
relaxation case as h = 0.01. The system size is there- 
fore scaled down to L^ x Ly x L^ = 64 x 64 x 64 to 
compensate the increased computational demand due to 
the small h explained above. The simulation results for 
e = 0.1 are shown in Fig.[5l which shows good agreement 
with the analytical solution of the relaxation function. 
The consistency between the input and calculated tem- 
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FIG. 3. (Color online) The fluid velocity field and its rota- 
tion, V X v, at f/r^ = 0.27. The values are normalized by the 
factor Ma* /\P\. The component normal to the figure plane 
is depicted. The compressibility factors are (a) e = 0.1 and 
(b) e — 1.5. The color scale (gray scale) represents negative 
(darker) to positive (lighter) rotation. The black circle repre- 
sents the particle. The direction of the initial particle velocity 
is right in these pictures. 
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FIG. 4. (Color online) The fluid velocity field and its diver- 
gence, V-i;, at t/Tv = 0.27. The values are normalized by the 
factor Ma* /\P\. The compressibility factors are (a) e = 0.1 
and (b) e — 1.5. The color scale (gray scale) represents nega- 
tive (darker) to positive (lighter) divergence. The black circle 
represents the particle. The direction of the initial particle 
velocity is right in these pictures. 



-is; 

en -1 ^ 

> 10 ^ 







10 



- £=0.1 

analytical solution (incompressible) 
analytical solution (£=0.1) 



J I I I I I I I I I I u 



10 



10' 



10" 
t/T„ 



10 



10 



FIG. 5. (Color online) The velocity autocorrelation function 
at e = 0.1 and ksT = 10~*. The bold solid line represents 
the analytic solution of the velocity relaxation function. The 
bold dashed double-dotted line shows the analytic solution for 
an incompressible fluid. 



Further improvenients on the treatments of fluctuations 
in the particle-fluid intcrfacial region will be discussed in 
the future. 



IV. CONCLUSION 

We extended the SPM to particle dispersions in com- 
pressible fluids. The validity of the method was con- 
firmed by calculating the velocity relaxation function of 
a single spherical particle in a compressible fluid. The 
effect of compressibility on the velocity relaxation was 
also observed, showing two-stage relaxation in a low- 
compressibility fluid and backtracking motion in a high- 
compressibility fluid. These particle motions were con- 
sidered by investigating the fluid density deviation. The 
propagation of the sound pulse around the particle is 
governed by the compressibility, and the influence of the 
sound disappears in a low-compressibility fluid but is 
maintained in a high-compressibility fluid. The effect 
of compressibility on the fluid velocity field was also ob- 
served, which was essentially understood to arise from 
changes in the time evolution of the source-sink flow com- 
ponent caused by the compressibility. 

A simulation of the motion of a single spherical particle 
in a fluctuating fluid was also performed. The calculated 
velocity autocorrelation function of the particle showed 
good agreement with the analytical solution of the re- 
laxation function, and the validity of the fluctuation- 
dissipation theorem was conflrmed. 
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peratures is also tested. We can evaluate the temperature 
from the average kinetic energy of a fluid or a particle, 
i.e., ksTf = A^(pv^)/3 or ksTp = M{V^)/3. We eval- 
uated these two temperatures as the ratio to the input 
temperature T. From the fluid motion, the temperature 
was evaluated as Tf/T ~ 1.06, while that of the parti- 
cle motion was Tp/T = 0.93. The overestimation of the 
fluid temperature Tf is simply due to the truncation er- 
ror in the time integration of the stochastic differential 
equation. On the other hand, the particle temperature 
Tp is slightly below the input temperature T. This dis- 
crepancy is considered to be due to the small numerical 
inaccuracy introduced in treating the momentum trans- 
fer through the particle-fluid interface using the SPM. 



APPENDIX: VELOCITY RELAXATION OF A 
PARTICLE IN A COMPRESSIBLE FLUID 

The equation of motion for a spherical particle under 
an external force E{t) in a compressible fluid can be ex- 
pressed by the linear response theory as 



dt 



dsCit - s)V{s) + E{t), (Al) 



where C(t) is the memory kernel of the friction force. The 
memory kernel is analytically expressed in frequency rep- 
resentation. Assuming a stick boundary condition on the 
surface of the sphere, the memory kernel is obtained from 



the linearized hvdrodynamic equations (Stokes approxi- 
mation) as [3,[a| 



C(^) 



with 



and 



((Oe'^^'di, 



iuj fi 






1/2 



(A2) 



— riax 
3 ' 

2x^(1 — iy) — (1 + x)y'^ — x^y^ 



X = a{-iLjpo/riy^^, y^auj/c, (A4) 



(A5) 



According to Eq. (jAip . the particle velocity is linearly 
dependent on the external force in frequency representa- 
tion: 



V{lu) = f{uj)E{uj), 
where the admittance r(w) is given by 

f(a;) = [-iujM + C{uj)]-\ 



(A6) 



(A7) 



In the case that an impulsive force is exerted on the 
sphere, the external force in frequency representation is 
given as a constant vector E{oj) = P, and the velocity 
relaxation function is given by 7(i) = Mr(i) according 
to Eq. (|22|) . 
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